rng('default'); 

% simulated state probabilities
prob_sim = mean(states(:)==(1:nS))';


% equity option IV
nKS    = 31;
KS_vec = linspace(-2,1,nKS)';
IV_KS  = NaN(nS,nK,nKS);
ATM    = IV(:,:,1);
for iks=1:nKS
    money = exp(KS_vec(iks)*ATM/sqrt(12)); % moneyness levels
    main_options
    IV_KS(:,:,iks) = BSIV;
end


% average across simulated states
EIV = NaN(nKS,1);
for iks=1:nKS
    IV_fit   = griddedInterpolant({1:nS,k},IV_KS(:,:,iks),'linear');
    EIV(iks) = nanmean(IV_fit(kron(ones(length(k_sim(:))/length(states(:)),1),states(:)),k_sim(:)));
end



% price asset put

[Cheb_nodes,~] = qnwcheb(nE,-1,1);
a           = -5;
b           = 5;
Eshocks     = (b-a)*(Cheb_nodes'+1)/2 + a;
temp        = pi*(b-a)/(2*nE);
weights     = temp.*sqrt(1-Cheb_nodes'.^2);
probE       = normpdf(Eshocks).*weights;


% ATM asset put
KS_value = 1;
lamAx = lamA - (1-tau_c)*(1-tau_d); 
g = exp(mu_XQ + sig_X.*Eshocks);
ATM_Asset_Put = NaN(nS,1);
ATM_Asset_IV = NaN(nS,1);
for is=1:nS % today's state
    ATM_Asset_Put(is) = (Q(is,:) * max(0,lamAx(is)*KS_value - lamAx*g(is,:)) * probE')/Rf(is,1);
    PV_div = (g(is,:)* probE')/Rf(is,1);
    ATM_Asset_IV(is) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_value*lamAx(is),log(Rf(is,1))*12,1/12,vol)-ATM_Asset_Put(is),0.2);
end


% IV skew in ATM asset IV units
nKS = 31;
MON_values = linspace(-2,1,nKS)';
KS_values = exp(MON_values * (ATM_Asset_IV'/sqrt(12)));
Asset_Put = NaN(nS,nKS);
Asset_IV = NaN(nS,nKS);
for is=1:nS % today's state
    PV_div = (g(is,:)* probE')/Rf(is,1);
    for iks=1:nKS
        Asset_Put(is,iks) = (Q(is,:) * max(0,lamAx(is)*KS_values(iks,is) - lamAx*g(is,:)) * probE')/Rf(is,1);
        Asset_IV(is,iks) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_values(iks,is)*lamAx(is),log(Rf(is,1))*12,1/12,vol)-Asset_Put(is,iks),0.2);
    end
end


% Asset-to-earnings ratio
% eps_it : Idiosyncratic shock, corr with eps_ct is rho (under P and Q)
[Cheb_nodes,~]  = qnwcheb(nE,-1,1);
a               = -5;
b               = 5;
Eshocks         = (b-a)*(Cheb_nodes'+1)/2 + a; % Heer and Maussner p.600
temp            = pi*(b-a)/(2*nE);
weights         = temp.*sqrt(1-Cheb_nodes'.^2);
probE           = normpdf(Eshocks).*weights;
probE           = probE';
Eshocks         = Eshocks';

X           = (exp(mu_XQ + sig_XS.*Eshocks')*probE) ./ Rf(:,1); % [nS,1]
lamA_noIdio = (1-tau_c)*(1-tau_d)*((eye(nS)-X.*Q)\ones(nS,1)); % [nS,1]



[Cheb_nodes,~] = qnwcheb(nE,-1,1);
a           = -5;
b           = 5;
Eshocks     = (b-a)*(Cheb_nodes'+1)/2 + a;
temp        = pi*(b-a)/(2*nE);
weights     = temp.*sqrt(1-Cheb_nodes'.^2);
probE       = normpdf(Eshocks).*weights;



% ATM asset put
KS_value = 1;
lamAx = lamA_noIdio - (1-tau_c)*(1-tau_d); 
g = exp(mu_XQ + sig_XS.*Eshocks);
ATM_Asset_Put_noIdio = NaN(nS,1);
ATM_Asset_IV_noIdio = NaN(nS,1);
for is=1:nS % today's state
    ATM_Asset_Put_noIdio(is) = (Q(is,:) * max(0,lamAx(is)*KS_value - lamAx*g(is,:)) * probE')/Rf(is,1);
    PV_div = (g(is,:)* probE')/Rf(is,1);
    ATM_Asset_IV_noIdio(is) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_value*lamAx(is),log(Rf(is,1))*12,1/12,vol)-ATM_Asset_Put_noIdio(is),0.2);
end


% IV skew in ATM asset IV units
nKS = 31;
MON_values = linspace(-2,1,nKS)';
KS_values = exp(MON_values * (ATM_Asset_IV_noIdio'/sqrt(12)));
Asset_Put_noIdio = NaN(nS,nKS);
Asset_IV_noIdio = NaN(nS,nKS);
for is=1:nS % today's state
    PV_div = (g(is,:)* probE')/Rf(is,1);
    for iks=1:nKS
        Asset_Put_noIdio(is,iks) = (Q(is,:) * max(0,lamAx(is)*KS_values(iks,is) - lamAx*g(is,:)) * probE')/Rf(is,1);
        Asset_IV_noIdio(is,iks) = fzero(@(vol)BS_put(lamAx(is)-PV_div,KS_values(iks,is)*lamAx(is),log(Rf(is,1))*12,1/12,vol)-Asset_Put_noIdio(is,iks),0.2);
    end
end


